/*
written for 
K.S. Gomez-Haibach and M.A. Gomez: Revised Centrality Measures Tell a Robust Story of Ion Conduction in Solids. 
The Journal of Physical Chemistry B 127, 9258 (2023).
using kinetic Monte Carlo scheme in
A.F. Voter: Introduction to the kinetic Monte Carlo method, in Radiation Effects in Solids, edited by K. E. Sickafus, E. A. Kotomin and B. P. Uberuaga (Springer Netherlands, Dordrecht, 2007)
Code from Gomez group at Mount Holyoke College - adapted for this tutorial
*/
// Headers
#include<stdio.h>
#include<stdlib.h>
#include <math.h>

//Numerical Recipies

#include "nrutil.h"
#include "ran3.cpp"

int main(){
	
	//Initialization
	// reading in data
	
	FILE *SitesFile;
	SitesFile =fopen("Energies.binding","r"); 
	if (SitesFile ==0){
		fprintf(stderr,"Can't open File Energies.binding");
		exit(1);
	}
	
	FILE *ConnectionsFile;
	ConnectionsFile=fopen("Energies.TS","r");
	if (ConnectionsFile==0){
        fprintf(stderr,"Can't open File Energies.TS");
        exit(1);
	}
	
	int numberofsites;
	fscanf(stdin, "%d\n", &numberofsites); //reads number of sites from input file
	fprintf(stdout,"%d\n",numberofsites);
	
	//Setting up the temperature and beta
	//1.3806503e-23 is Boltzmann's constant in units of J/K
	//6.24150974e18 is the conversion from J to eV
	double beta, temperature, kb;
	kb=1.3806503e-23*6.24150974e18;
	fprintf(stdout,"What temperature would you like to run at? (K)\n");
	fscanf(stdin,"%lf\n",&temperature);
	fprintf(stdout,"%lf Kelvin\n",temperature);
	beta=1/(kb*temperature);
	fprintf(stdout,"beta is %lf in units of 1/eV\n",beta);
	
	double **Bpositions, *siteProbability, *Energy;
	Bpositions = dmatrix(1, numberofsites, 1, 3);
	Energy=dvector(1,numberofsites);
	siteProbability = dvector(1, numberofsites);
	
	int i, idum, j, jdum;
	
	for(i=1; i<=numberofsites; i++){
		fscanf(SitesFile, "%d %lf %lf %lf %lf\n", &idum, &Bpositions[i][1], &Bpositions[i][2], &Bpositions[i][3], &Energy[i]);
		//fprintf(stdout, "%d %lf %lf %lf %lf\n", idum, Bpositions[i][1], Bpositions[i][2], Bpositions[i][3], Energy[i]);
	}
	
	//Finding the lowest energy site
	double minEnergy;  minEnergy=0.0;//large number as all minima have negative energies
	int imin;
	for (i=1;i<=numberofsites;i++) {
			if (minEnergy>Energy[i]) {
				minEnergy=Energy[i];
				imin=i;
			}
	}
	fprintf(stdout,"Mininum Energy is %lf and occurs for binding site %d\n",minEnergy,imin);
	
	double boltznorm; boltznorm=0.0;
	//Boltmann normalization factor
	for (i=1;i<=numberofsites;i++) boltznorm=boltznorm+exp(-beta*(Energy[i]-minEnergy)); 
	// if we had frequencies, each term would need to be divided by the product of frequencies at the ith site
	//fprintf(stdout,"boltz norm %lf\n",boltznorm);
	
	//fprintf(stdout,"Site Position Energy Probability\n");
	for (i=1;i<=numberofsites;i++) {
		siteProbability[i]=exp(-beta*(Energy[i]-minEnergy))/boltznorm;// probability would also need to be divided by the product of frequencies at each state if we had frequencies
		//fprintf(stdout,"%d %e\n",i,siteProbability[i]);
	}
	double TSEnergy;
	
	double **Pij, **sumOfPij, **RateConstant, *ktot, *sumOfSiteProbabilities, **TS, *sumTimeMoveType,**MoveTimes;
	int **Move, *countMoveType;
	
	Pij = dmatrix(1, numberofsites, 1, numberofsites);
	sumOfPij = dmatrix(1, numberofsites, 1, numberofsites);
	RateConstant = dmatrix(1, numberofsites, 1, numberofsites);
	MoveTimes = dmatrix(1, numberofsites, 1, numberofsites);
	TS = dmatrix(1,numberofsites,1,numberofsites);
	Move = imatrix(1,numberofsites,1,numberofsites);
	countMoveType=ivector(1,3);  countMoveType[1]=0; countMoveType[2]=0; countMoveType[3]=0;
	sumTimeMoveType=dvector(1,3);  sumTimeMoveType[1]=0.0; sumTimeMoveType[2]=0.0; sumTimeMoveType[3]=0.0;
	ktot = dvector(1,numberofsites);
	sumOfSiteProbabilities = dvector(1,numberofsites);
	
	for (i=1; i<=numberofsites; i++){
		sumOfSiteProbabilities[i]=0.0;
		for (j=1;j<=i;j++) 
			sumOfSiteProbabilities[i]=sumOfSiteProbabilities[i]+siteProbability[j];
		for(j=1; j<=numberofsites; j++){
			Pij[i][j]=0.0;
			sumOfPij[i][j]=0.0;
			RateConstant[i][j]=0.0;
			TS[i][j]=0.0;
			MoveTimes[i][j]=0.0;
			Move[i][j]=0;
		}
		ktot[i]=0.0;
	}
	
	int numOfProbabilities,imove; // Number of non-zero probabilities.
	double xdum,ydum,zdum;
	
	numOfProbabilities=216; 
	for(i=1; i<=numOfProbabilities; i++){
		fscanf(ConnectionsFile, "%d %d %d %lf \n",&idum,&jdum,&imove,&TSEnergy); //reads transition state energies 
		RateConstant[idum][jdum]=exp(-beta*(TSEnergy-Energy[idum]));//would need to multiplied by a prefactor (ratio of frequencies at minimum to those at TS if we had frequencies
        	RateConstant[jdum][idum]=exp(-beta*(TSEnergy-Energy[jdum]));
		TS[idum][jdum]=TSEnergy; TS[jdum][idum]=TSEnergy;
		Move[idum][jdum]=imove; Move[jdum][idum]=imove;
		//fprintf(stdout, "%d %d %d %lf %lf %lf %lf %d %d\n",i,idum,jdum,xdum,ydum,zdum,TSEnergy,Move[idum][jdum],Move[jdum][idum]); //reads transition state energies 
	}

	//calculating Pij
	double sum;
	fprintf(stdout,"site probability, rateconstant\n");
	for(i=1;i<=numberofsites;i++){
		sum=0.0;
		for (j=1;j<=numberofsites;j++) sum=sum+RateConstant[i][j];
		for(j=1;j<=numberofsites;j++) {
			if (sum>0.0) {Pij[i][j]=RateConstant[i][j]/sum;} else {Pij[i][j]=0.0;}
			//if (Pij[i][j]>0) fprintf(stdout,"%d %d %e %e\n",i,j,Pij[i][j],RateConstant[i][j]);
		}
	}
	
	//calculate sumOfRateConstants before we start KMC
	int k;
	//fprintf(stdout, "i ktot sumPij starting at site i\n");
	for (i=1; i<=numberofsites; i++){ //position number for sum
		for(j=1; j<=numberofsites; j++){ //number of terms added in the sum
			for (k=1; k<=j; k++) sumOfPij[i][j]+=Pij[i][k];
			ktot[i]+=RateConstant[i][j];
		}
		//fprintf(stdout, "%d %lf %lf\n", i, ktot[i],sumOfPij[i][numberofsites]);
	}
	
	//setting up first site

	int iseed;	

	fscanf(stdin, "%d\n", &iseed);
  fprintf(stdout, "The seed is %d\n",iseed);
  fprintf(stdout,"sample random number %lf\n",ran3(iseed));

  double randomNumber, avgRandom, avg2Random;
  avgRandom=0.0; avg2Random=0.0;
  int npts=1000000;
  for (i=1;i<=npts;i++) {randomNumber=ran3(iseed); avgRandom=avgRandom+randomNumber; avg2Random=avg2Random+randomNumber*randomNumber;}
  fprintf(stdout,"Testing random number generator (comparisons via integration in parenthesis)\n");
  fprintf(stdout,"<x>=%lf (1/2) <x2>=%lf (1/3) npts=%d\n",avgRandom/double(npts),avg2Random/double(npts),npts);

	double targetTime;
	fscanf(stdin,"%lf\n",&targetTime);
	fprintf(stdout,"target time=%lf\n",targetTime);

	int startingSite,currentSite, prevSite, maxI, maxF;
	fscanf(stdin,"%d\n",&startingSite);
	
	double time,timePerCrossing;
	time=0.0;
	fprintf(stdout,"Starting\nTime Site Siteprobability\n%lf %d %lf\n",time, startingSite,siteProbability[startingSite]);
  fprintf(stdout,"Each tractory is stopped when it has reached the other side of the box\nHere are the final times, sites, maxBarrier and the location of it\n");

	double timestep;
	double randNum;
	
	double *Distribution;
	Distribution=dvector(1,numberofsites);
	
	double xdisp, ydisp, zdisp, dx, dy, dz, MaxBarrier,currentBarrier,avgMaxBarrier,avgMaxBarrierTime;
	int flag,numberCrossings,ncross;
  //change the number of crossings here
	numberCrossings=1000;
	avgMaxBarrier=0.0;
	avgMaxBarrierTime=0.0;

//start of numberCrossings trajectories
fprintf(stdout,"\nIndividual trajectory information\nCrossingTime FinalSite MaxBarrierOnPath StartingAt EndingAt ViaMoveType(1 is R,2 is T,3 is I)\n");
for (ncross=1;ncross<=numberCrossings;ncross++) {

  //variables to keep track of how far a trajectory has moved
	xdisp = 0.0;
	ydisp = 0.0;
	zdisp = 0.0;
	flag = 0;
	MaxBarrier=0.0;
	timePerCrossing=0.0;
	currentSite=startingSite;  
	prevSite = currentSite;
 
  //loop to evolve trajectory until either crossingoccurs or target time runs out.  The flag changes from 0 to 1 when a crossing occurs
	while (timePerCrossing < targetTime && flag == 0){
 
		dx=Bpositions[currentSite][1] - Bpositions[prevSite][1];
		if (dx<-0.5) dx=dx+1.0;
 	  if (dx>0.5) dx=dx-1.0;
		xdisp = xdisp + dx;
		
		dy=Bpositions[currentSite][2] - Bpositions[prevSite][2];
		if (dy<-0.5) dy=dy+1.0;
 	  if (dy>0.5) dy=dy-1.0;
		ydisp = ydisp + dy;
		
		dz=Bpositions[currentSite][3]-Bpositions[prevSite][3];
		if (dz<-0.5) dz=dz+1.0;
 	  if (dz>0.5) dz=dz-1.0;
		zdisp = zdisp + dz;
		
    //checking if box is crossed
		if (xdisp > 0.95 || ydisp > 0.95 || zdisp > 0.95 || xdisp < -0.95 || ydisp < -0.95 || zdisp < -0.95){
   
			flag = 1; //indicates crossing or having traversed more than 95 percent of the box in any direction
      //fprintf(stdout, "time=%lf currentSite=%d xdisp=%lf ydisp=%lf zdisp=%lf flag=%d maxBarrier=%lf from %d to %d\n", time,currentSite,xdisp, ydisp, zdisp, flag, MaxBarrier,maxI,maxF);
		} else {
	
			//KMC steps continue if no crossing has occurred
      
			prevSite = currentSite;
      randNum=ran3(iseed); //random number between 0 and 1 is generated
      while (randNum==0 || randNum==1) randNum=ran3(iseed); //for the unlikely event that the random number is 0 or 1

      timestep=-log(randNum)/ktot[currentSite];  //time step calculated for move
      timePerCrossing=timePerCrossing+timestep; //updates time to keep track of how long it takes to cross in this trajectory
      Distribution[currentSite]=Distribution[currentSite]+timestep; //adds time spent at currentSite before move to monitor time at site

      randNum=ran3(iseed);
      while (randNum==0 || randNum==1) randNum=ran3(iseed); //for the unlikely event that the random number is 0 or 1

      j=1;
      while (randNum > sumOfPij[currentSite][j]) j++;  //finds the site probability section that contains the random number
      currentSite = j; //sets the new position

      //checking to see if the new current barrier is larger than the old max to eventually store the max barrier for this long
      //range path (accross box) and the transition involved
			currentBarrier=TS[prevSite][currentSite]-Energy[prevSite];
			if (currentBarrier>MaxBarrier) {MaxBarrier=currentBarrier; maxI=prevSite; maxF=currentSite;}	
      //fprintf(stdout, "time=%lf currentSite=%d xdisp=%lf ydisp=%lf zdisp=%lf flag=%d currentBarrierr=%lf\n", time,currentSite,xdisp, ydisp, zdisp, flag, currentBarrier);
		}
	}
 
  //summing the time for the crossings, the max barriers, the max barrier times the time ... to get ready for averages
	time=time+timePerCrossing;
	avgMaxBarrier=avgMaxBarrier+MaxBarrier;
	avgMaxBarrierTime=avgMaxBarrierTime+MaxBarrier*timePerCrossing;
 
  //printing individual trajectory information
  fprintf(stdout, "%lf %d maxBarrier=%lf from %d to %d via %d\n", timePerCrossing,currentSite,MaxBarrier,maxI,maxF,Move[maxI][maxF]);
	countMoveType[Move[maxI][maxF]]=countMoveType[Move[maxI][maxF]]+1;
	sumTimeMoveType[Move[maxI][maxF]]=sumTimeMoveType[Move[maxI][maxF]]+timePerCrossing;
	MoveTimes[maxI][maxF]=MoveTimes[maxI][maxF]+timePerCrossing;
	if (flag!=1) fprintf(stdout,"Increase maximum time to see a move accross the simulation box.\n");
 
}

  //average calculations and printing
  fprintf(stdout,"\n\nAverages over trajectories that cross simulation box\n");
	fprintf(stdout,"avgMaxBarrier(over number of crossings, over time)=(%lf,%lf) avgTimePerCrossing=%lf\n",avgMaxBarrier/double(numberCrossings),avgMaxBarrierTime/time,time/double(numberCrossings));
  fprintf(stdout,"  fractionOfCrossings fractionOfTime\n");
	fprintf(stdout,"R %lf,%lf \nT %lf,%lf \nI %lf,%lf\n",double(countMoveType[1])/double(numberCrossings),sumTimeMoveType[1]/time,double(countMoveType[2])/double(numberCrossings), sumTimeMoveType[2]/time,double(countMoveType[3])/double(numberCrossings),sumTimeMoveType[3]/time);

	fprintf(stdout,"\n\nIndividual Moves\nfrom to MoveType timeFractionOfMove\nIf you want these printed, uncomment lines.\n");
	/*for (i=1;i<=numberofsites;i++) {
		for (j=1;j<=numberofsites;j++) {
			if (MoveTimes[i][j]>0.) fprintf(stdout,"%d %d %d %lf\n",i,j,Move[i][j],MoveTimes[i][j]/time);
		}
	}
 */
	
	fprintf(stdout,"\nCheck that time fraction and site probability are about the same to insure that time is long enough and/or that there are enough trajectories\nx y z j HNumber TimeFraction=fractionOfTheTimeSpentAtThisLocation SiteProbability\nIf you want these printed, uncomment lines.\n");
 
 /*
	for (i=1;i<=numberofsites;i++) {
		if (Distribution[i]>0) {
			fprintf(stdout, "%lf %lf %lf H%d TimeFraction=%e siteProbability=%e\n",Bpositions[i][1],Bpositions[i][2],Bpositions[i][3],i,Distribution[i]/time,siteProbability[i]);
		}
	}*/
 
}
